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Coexistence versus extinction in the stochastic cyclic Lotka-Volterra model 
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Cyclic dominance of species has been identified as a potential mechanism to maintain biodiversity, 
see e.g. B. Kerr, M. A. Riley, M. W. Feldman and B. J. M. Bohannan [Nature 418, 171 (2002)] 
and B. Kirkup and M. A. Riley [Nature 428, 412 (2004)]. Through analytical methods supported 
by numerical simulations, we address this issue by studying the properties of a paradigmatic non- 
spatial three-species stochastic system, namely the 'rock-paper-scissors' or cyclic Lotka-Volterra 
model. While the deterministic approach (rate equations) predicts the coexistence of the species 
resulting in regular (yet neutrally stable) oscillations of the population densities, we demonstrate 
that fluctuations arising in the system with a finite number of agents drastically alter this picture 
and are responsible for extinction: After long enough time, two of the three species die out. As main 
findings we provide analytic estimates and numerical computation of the extinction probability at a 
given time. We also discuss the implications of our results for a broad class of competing population 
systems. 

PACS numbers: 05.40.-a,87.23.Cc,02.50.Ey,05.10.Gg 



I. INTRODUCTION 



Understanding biodiversity and coevolution is a central 
challenge in modern evolutionary and theoretical biology 
01 EH- In this context, for some decades much effort has 
been devoted to mathematically model dynamics of com- 
peting populations through nonlinear, yet deterministic, 
set of rate equations like the equations devised by Lotka 
and Volterra 0, or many of their variants 0, 0, Q- 
This heuristic approach is often termed as population- 
level description. As a common feature, these determin- 
istic models fail to account for stochastic effects (like fluc- 
tuations and spatial correlations). However, to gain some 
more realistic and fundamental understanding on generic 
features of population dynamics and mechanisms leading 
to biodiversity, it is highly desirable to include internal 
stochastic noise in the description of agents' kinetics by 
going beyond the classical deterministic picture. One of 
the main reasons is to account for discrete degrees of free- 
dom and finite-size fluctuations 0,0- In fact, the deter- 
ministic rate equations always (tacitly) assume the pres- 
ence of infinitely many interacting agents, while in real 
systems there is a large, yet finite, number of individuals 
(recently, this issue has been addressed in Refs. |^,@). 
As a consequence, the dynamics is intrinsically stochas- 
tic and the unavoidable finite-size fluctuations may have 
drastic effects and even completely invalidate the deter- 
ministic predictions. 

Interestingly, both in vitro |lOj and in vivo fllT ] ex- 
periments have recently been devoted to experimentally 
probe the influence of stochasticity on biodiversity: The 
authors of Refs. [ToL ITlT j have investigated the mecha- 
nism necessary to ensure coexistence in a community of 
three populations of Escherichia coli and have numeri- 
cally modelled the dynamics of their experiments by the 
so-called 'rock-paper-scissors' model, well-known in the 



field of game theory [H[l3. This is a three-s pec ies cyclic 
generalization of the Lotka-Volterra model |14l fl5L fl6| . 

As a result, the authors of Ref. reported that in a 
well-mixed (non-spatial) environment (i.e. when the ex- 
periments were carried out in a flask) two species got 
extinct after some finite time, while coexistence of the 
populations was never observed. Motivated by these ex- 
perimental results, in this work we theoretically study the 
stochastic version of the cyclic Lotka- Volterra model and 
investigate in detail the effects of finite-size fluctuations 
on possible population extinction/coexistence. 

For our inv estig ation, as suggested by the flask exper- 
iment of Ref. [l(| , the stochastic dynamics of the cyclic 
Lotka-Volterra model is formulated in the natural lan- 
guage of urn models 01 and by adopting the so-called 
individual-based description [lq. In the latter, the ex- 
plicit rules governing the interaction of a finite number 
of individuals with each other are embodied in a mas- 
ter equation. The fluctuations are then specifically ac- 
counted for by an appropriate Fokker-Planck equation 
derived from the master equation via a so-called van 
Kampen expansion |l9l] . This program allows us to quan- 
titatively study the deviations of the stochastic dynam- 
ics of the cyclic Lotka-Volterra model with respect to the 
rate equation predictions and to address the question of 
the extinction probability, the computation of which is 
the main result of this work. From a more general per- 
spective, we think that our findings have a broad rele- 
vance, both theoretical and practical, as they shed fur- 
ther light on how stochastic noise can dramatically affect 
the properties of the numerous nonlinear systems whose 
deterministic description, like in the case of the cyclic 
Lotka-Volterra model, predicts the existence of neutrally 
stable solutions, i.e. cycles in the phase portrait [20| . 

This paper is organized as follows: The cyclic Lotka- 
Volterra model is introduced in the next section and its 
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FIG. 1: Illustration of cyclic dominance of three states A, 
B, and C. The latter m ay c orrespond to the strategies in a 
rock-paper-scissors game or to different bacterial species 
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deterministic rate equation treatment is presented. In 
Section lTTTl we develop a quantitative analytical approach 
that accounts for stochasticity, a Fokker-Planck equation 
is derived from the underlying master equation within 
a Van Kampen expansion. This allows us to compute 
the variances of the agents' densities. We also study the 
time-dependence properties of the system by carrying out 
a Fourier analysis from a set of Langevin equations. Sec- 
tion ^^D] is devoted to the computation of the probabil- 
ity of having extinction of two species at a given time, 
which constitutes the main issue of this work. In the fi- 
nal section, we summarize our findings and present our 
conclusions. 



II. THE CYCLIC POPULATION MODEL: 
INTRODUCTION AND ANALYSIS OF THE 
RATE EQUATIONS 

A. The model 

The cyclic Lotka Volterra model under consideration 
here is a system where three states A, B, and C cycli- 
cally dominate each other: A invades B, B outperforms 
C, and C in turn dominates over A, schematically drawn 
in Fig. 2] These three states A, B, and C allow for var- 
ious interpretations, ranging from strategies in the rock- 
paper-scissors game 01 over tree, fire, ash in forest fire 
models [2lJ or chemical reactions [2^ to different bac- 
terial species 0, In the latter case, a population 
of poison producing bacteria was brought together with 
another being resistant to the poison and a third which 
is not resistant. As the production of poison as well as 
the resistance against it have some cost, these species 
show a cyclic dominance: the poison-producing one in- 
vades the non-resistant, which in turn reproduces faster 
than the resistant one, and the latter finally dominates 
the poison-producing. In a well-mixed environment, like 
a flask in the experiments [ToL ITlj , eventually only one 
species survives. The populations are large but finite, 
and the dynamics of reproduction and killing events may 
to a good approximation be viewed as stochastic. 

Motivated from these biological experiments, we in- 
troduce a stochastic version of the cyclic Lotka- Volterra 
model. Consider a population of size ./V which is well 
mixed, i.e. in the absence of spatial structure. The 





FIG. 2: (color online) The urn model describes stochastic 
evolution of well-mixed finite populations. We show three 
species as yellow or light gray (A), red or medium gray (B), 
and blue or dark gray (C). At each time step, two random 
individuals are chosen (indicated by arrows in the left picture) 
and react (right picture). 



stochastic dynamics used to describe its evolution, illus- 
trated in Fig. [2 is referred to as "urn model" [17| and 
closely related to the Moran process [2^, [3] • At every 
time step, two randomly chosen individuals are selected, 
which may at certain probability react according to the 
following scheme: 
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with reaction rates fc^, ks, and kc- We observe the 
cyclic dominance of the three species. Also, the total 
number N of individuals is conserved by this dynamics; 
this will of course play a role in our further analysis. 

We now proceed with the analysis of the deterministic 
version of the system QJ. This will prove insightful for 
building a stochastic description of the model, which is 
the scope of Sec. IIIII 



B. Deterministic approach 

The deterministic rate equations describe the time evo- 
lution of the densities a(t), b(t), and c(t) for the species 
A, B, and C; they read 



a = a(kcb — ksc) 
b — b(kAC — kca) 
c = c(ksa — k^b), 



(2) 



where the dot stands for the time derivative. These equa- 
tions describe a well-mixed system, without any spatial 
correlations, as naturally implemented in urn models 0] 
or, equivalently, infinite dimensional lattice systems or 
complete graphs. In the following, the Eqs. are dis- 
cussed and, from their properties, we gain intuition on 
the effects of stochasticity. 
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Already from the basic reactions Q we have noticed 
that the total number of individuals is conserved, which is 
a property correctly reproduced by the the rate equations 
(J2J) ■ Setting the total density, meaning the sum of the 
densities a, b, and c, to unity, we obtain 



a{t) + b(t) + c(t) = 1 



(3) 



for all times t. Only two out of the three densities are 
thus independent, we may view the time evolution of the 
densities in a two-dimensional phase space. 

Eqs. (0) together with J2J admit three trivial (absorb- 
ing) fixed points: (a*,b*,c*) = (1,0,0); (a^b^ci,) = 
(0,1,0); and (a%,b%,c%) = (0,0,1). They denote states 
where only one of the three species survived, the other 
ones died out. In addition, the rate equations also 
predict the existence of a fixed point (a*,6*,c*) which 
corresponds to a reactive steady state, associated with 
the coexistence of all three species: 
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(4) 



kA + ks + kc 

To determine the nature of this fixed point, we ob- 
serve that another constant of motion exists for the rate 
equations iffil. namely the quantity 

a{t) kA b(t) kB c(t) kc = a{0) kA b(0) kB c(0) kc (5) 

does not evolve in time. In contrast to the total density 
iJSJ , this constant of motion is only conserved by the rate 
equations but does not stem from the reaction scheme 
JTJ. Hence, when considering the stochastic version of 
the cyclic model, the total density remains constant but 
the expression J5J will no longer be a conserved quantity. 
The above fixed point Q and constant of motion JSJ have 
been derived and discussed also within the framework of 
game theory, see e.g. Ref. [T^ . 

In Fig. |2| we depict the ternary phase space [2f| for 
the densities a, &, and c: The solutions of the rate equa- 
tions (|2J) are shown for different initial conditions and a 
given set of rates kA, fcs, and kc- As the rate equa- 
tions are nonlinear in the densities, we cannot solve 
them analytically, but use numerical methods. Due to 
the constant of motion, the solutions yield cycles around 
the reactive fixed point (thus corresponding to case 3 in 
Durrett and Levin's classification |2(j ) . The three trivial 
steady states, corresponding to saddle points within the 
linear analysis, are the edges of the simplex. The reac- 
tive stationary state, as well as the cycles, are neutrally 
stable, stemming from the existence of the constant of 
motion |JSj. Especially, the reactive fixed point is a cen- 
ter fixed point. The boundary of the simplex denotes 
states where at least one of the three species died out; 
as cyclic dominance is lost, states that have reached this 
boundary will evolve towards one of the edges, making 
the boundary absorbing. 




FIG. 3: (color online) The three-species simplex for reaction 
rates kA = 1, fee = 1-5, kc = 2. The rate equations predict 
cycles, which are shown in black. Their linearization around 
the reactive fixed point is solved in proper coordinates j/a, Vb 
(blue or dark gray). The red (or light gray) erratic flow de- 
notes a single trajectory in a finite system (JV = 200), ob- 
tained from stochastic simulations. It spirals out from the 
reactive fixed point, eventually reaching an absorbing state. 



1. Linear regime 

The nonlincarity of Eqs. @ induces substantial diffi- 
culties in the analytical treatment. However, much can 
already be inferred from the linearization around the re- 
active fixed point which we will consider in the fol- 
lowing. We therefore introduce the deviations from the 
reactive fixed point, denoted as, xa, xb, xc- 



xa = a — a 
xb = b — b* 
xc = c — c* 



(6) 



Using conservation of the total density © , we can elim- 
inate xc-, and the remaining linearized equations @ 
may be put into the form x = Ax, with the vector 
x = (xa,xb) and the matrix 



A = 



1 



k A + k B + k c 



k A k B k A (k B + kc) 

-k B (k A + kc) -k A k B 



(7) 

The reactive fixed point is associated to the eigenvalues 
A± = iio^o of A, where luq is given by 



UJ 



kAksK 



kA + ks + kc 



(8) 



oscillations with this frequency arise in its vicinity. In 
proper coordinates y ~ <Sx, these oscillations are har- 
monic, being the solution of y = Ay, with A = 

SAS^ 1 = f ° U n V For illustration, we have included 

these coordinates, in which the solutions take the form 
of circles around the origin, in Fig. [31 The linear trans- 
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FIG. 4: (color online) The deterministic time-evolution of 
the density a for small and large amplitudes. The prediction 
(lilt , shown in black, is compared to the numerical solution 
(red or gray) of the rate equations (j2J. For small amplitudes 
(a(0) = a* + 0.03, &(0) = b*), both coincide. However, for 
large amplitudes (a(0) = a* + 0.2, 6(0) = 6*), they consider- 
ably differ both in amplitude and frequency. We used reaction 
rates h,A = 1, fcs = 1.5, kc = 2. 



formation x 



S 



y is given by 
_ V3 (kA+kc 



2 V 



1 



(9) 



The equations y = Ay are easily solved: 

y A (t) = y A (0) cos(^oi) + Vb{0) sin(u> t) 
VB(t) = y B {0) cos(uj t) - y A (0) sin(w f) 



(10) 



Eventually, we obtain the solutions for the linearized rate 
equations: 



o(t) 



[o(0) 



i ^[o(0) 



o*] cos(woi) 

+ fcg + fcc 

k B k c 



6(0) 



sin(wot) 
(11) 



where b(t) and c(i) follow by cyclic permutations. 

To establish the validity of the linear analysis (|7|)- (|llfl . 
we compare the (numerical) solution of the rate equa- 
tions with the linear approximations l|llf> . As shown 
in Fig. H when o(0) - a* « 0.03, 6(0) - 6* = 0, the 
agreement between the nonlinear rate equations (0 and 
the linear approximation (|llf> is excellent, both curves 
almost coincide. On the other hand, the nonlinear terms 
appearing in Eqs. @ become important already when 
a(0) — a* « 0.1 and are responsible for significant dis- 
crepancies both in the amplitudes and frequency from 
the predictions of Eq. (|llfl . 



2. Radius TZ 

We now aim at introducing a measure of distance TZ 
to the reactive fixed point within the phase portrait. In 
the next section, this quantity will help quantify effects of 
stochasticity. As it was recently proved useful in a related 
context ja], we aim at taking the structure of cycles 



predicted by the deterministic equations, see Fig. |3J into 
account by requiring that the distance should not change 
on a given cycle. Motivated from the constant of motion 
(151). we introduce 



TZ = Af\/a* kA b* kB c* kc -a kA b ks c kc 
with the normalization factor (see below) 



k B {k A + kc) 



2 (k A + k B + k c ) 3 a* kA b* kB c* k c 



(12) 



(13) 



Being conserved by the Eqs. J2J), TZ remains constant on 
every deterministic cycle. As it vanishes at the reac- 
tive fixed point and monotonically grows when departing 
from it, TZ yields a measure of the distance to the latter. 

Expanding the radius TZ in small deviations x from the 
reactive fixed point results in 



A" 



ft =^ r (k A + k B + kcYa* KA b 



h> 2A+ (i- B + i> 



kc 



x A x B 



+ o(x 2 ) 



(14) 



In the variables y = iSx, with our choice for J\f, it sim- 
plifies to 



7^- 



y A + y 2 B + o(y 2 ) 



(15) 



corresponding to the radius of the deterministic circles, 
which emerge in the variables y. 



3. Effects of stochasticity: qualitative discussion 

The fact that the number N of particles is finite in- 
duces fluctuations that are not accounted for by Eqs. J5J). 
In the following, our goal is to understand the impor- 
tance of fluctuations and their effects on the determin- 
istic picture (J2J. We show that, due to the neutrally 
stable character of the deterministic cycles, fluctuations 
have drastic consequences. Intuitively, we expect that in 
the presence of stochasticity, each trajectory performs a 
random walk in the phase portrait, interpolating between 
the deterministic cycles (as will be revealed by consider- 
ing TZ), eventually reaching the boundary of the phase 
space. There, the cyclic dominance is completely lost, 
as one of the species gets extinct. Of the two remaining 
ones, one species is defeating the other, such that the lat- 
ter soon gets extinct as well, leaving the other one as the 
only survivor (this corresponds to one of the trivial fixed 
points, the edges of the ternary phase space). We thus 
observe the boundary to be absorbing, and presume the 
system to always end up in one of the absorbing states. 
A first indication of the actual emergence of this scenario 
can be inferred from the stochastic trajectory shown in 
Fig. 
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III. FLUCTUATIONS IN FINITE 
POPULATIONS: A STOCHASTIC APPROACH 

In this section, we set up a stochastic description for 
the cyclic Lotka-Voltcrra model in the urn formulation 
with finite number N of individuals. Starting from the 
master equation of the stochastic process, we obtain a 
Fokker-Planck equation for the time evolution of the 
probability P(a, b, c, t) of finding the system in state 
a, 6, c at time t. It allows us to gain a detailed under- 
standing of the stochastic system. In particular, wc will 
find that, as anticipated at the end of the last section, 
after long enough time, the system reaches one of the ab- 
sorbing states. Our main result is the time-dependence 
of the extinction probability, being the probability that, 
starting at a situation corresponding to the reactive fixed 
point, after time t two of the three species have died 
out. It is obtained through mapping onto a known first- 
passage problem. Wc compare our analytical findings to 
results from stochastic simulations. For the sake of clar- 
ity and without loss of generality, throughout this sec- 
tion, the case of equal reaction rates kj± — kg = kc = 1 
is considered. Details on the unequal rates situation are 
relegated to Appendix IaI 

A. Stochastic Simulations 

We carried out extensive stochastic simulations to sup- 
port and corroborate our analytical results. An efficient 
simulation method originally due to Gillespie [27J was 
implemented for the reactions ffl. Time and type of 
the next reaction taking place are determined by random 
numbers, using the Poisson nature of the individual reac- 
tions. For the extinction probability, to unravel the uni- 
versal time-scaling, system sizes ranging from ./V = 100 
to TV = 1000 were considered, with sample averages over 
10000 realizations. 



B. From the master to Fokker-Planck equation 

Let us start with the master equation of the processes 
||TJ. We derive it in the variables xa,%b>%C> which were 
introduced in (JfJJ) as the deviations of the densities from 
the reactive fixed point. Using the conservation of the 
total density, xa + xb + %c = 0- xc is eliminated and 
x = (xa,xb) kept as independent variables. The Master 
equation for the time-evolution of the probability P(x, t) 
of finding the system in state x at time t thus reads 

d t P(x, t) = X^{ p ( x + <5x > *) W ( X + Sx -> x ) 

<5x 

-P(x,i)W(x^x + <Sx)} , (16) 

where W(x — > x + <5x) denotes the transition probability 
from state x to the state x+<5x within one time step; sum- 
mation extends over all possible changes (5x. Wc choose 



the unit of time so that, on average, every individual re- 
acts once per time step. 

According to the Kramers-Moyal expansion of the Mas- 
ter equation to second order in 5x, this results in the 
Fokker-Planck equation 

d t P(x,t) = -5 i [a i (x)P(x,t)] + ifta,[S y -(x)P(x,t)] . 

(17) 

Here, the indices i,j stand for A and B; in the above 
equation, the summation convention implies summation 
over them. The quantities cdj and Bij arc, according to 
the Kramcrs-Moyal expansion: 

a; (x) = J~] faj>V(x — > x + 5x) 

<5x 

B i3 {x) =J25x i 5x J W(x^x + Sx) . (18) 

,5x 

Note that B is symmetric. For the sake of clarity, wc 
outline the calculation of 04 (x): The relevant changes 
Sxa in the density a result from the basic reactions 
||TJ, they are 6xa = 1/iV in the first reaction, Sxa = 
in the second and Sxa = —1/N in the third. The 
corresponding rates read W = Nab for the first reaction 
(the prefactor of N enters due to our choice of time 
scale, where N reactions occur in one unit of time), and 
VV = Nac for the third, resulting in Q!a(x) = ab — ac. 
The other quantities are calculated analogously, yielding 

a,! (x) = ab — ac 

ag(x) = be — ab 
&aa( x ) = N~ 1 (ab + ac) (19) 
Bab{x) = B B a{x) = -N^ab 
S BB (x) = N-^bc + ab) 

with a = a* + xa, b = b* + xb, c = c* — xa — xb- 

Van Kampen's linear noise approximation |l9j further 
simplifies these quantities. In this approach, the values 
a,(x) are expanded around x = to the first order. As 
they vanish at the reactive fixed point, we obtain: 

a-j(x) = AijXj, (20) 

where the matrix elements Aij are given by 

•Aii = ^ • (21) 

axj x=o 

The matrix A, already given in Q , embodies the deter- 
ministic evolution, while the stochastic noise is encoded 
in B. To take the fluctuations into account within the 
Van Kampen expansion, one approximates B by its val- 
ues at the reactive fixed point. Hence, we find: 
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and the corresponding Fokker-Planck equation reads 
$P(x,t) = -d i [A ij x j P{xL,t)]+^B ij d i d j P(x,t) . (23) 

For further convenience, we now bring Eq. (|23[) into a 
more suitable form by exploiting the polar symmetry 
unveiled by the variables y. As for the linearization 
of Eqs. (J2J), it is useful to rely on the linear mapping 
x — * y = <Sx, with A — > A = <SAS _1 . Interestingly, it 
turns out that this transformation diagonalizes B. One 
indeed finds B — > B = SBS T , with 



B 



J_ (l 

6JV lo 1 



(24) 



In the y variables, the Fokker-Planck equation i|23|) takes 
the simpler form 

dtP{y, t)=- Lo Q [y A dy B - y B d yA ]P(y, t) 

1 M A +9L}P(y,t) ■ (25) 



127V 



To capture the structure of circles predicted by the de- 
terministic approach, we introduce polar coordinates: 



in the following. Within our approximations around the 
reactive fixed point, the probability distribution P(r, t) 
is thus the same as for a system performing a two- 
dimensional random walk in the y variables. Intuitively, 
such behavior is expected and its origin lies in the neu- 
trally stable character of the cycles of the deterministic 
solution of the rate equations J5J). The cycling around 
the reactive fixed point does not yield additional effects 
when considering spherically symmetric probability dis- 
tributions. Furthermore, due to the existence of the con- 
stant of motion JSJ), the neutral stability does not only 
hold in the vicinity of the reactive fixed point, but in the 
whole ternary phase space. We thus expect the distance 
TZ from the reactive fixed point, a random variable, to 
obey a diffusion-like equation in the whole phase space 
(and not only around the reactive fixed point). Actually, 
qualitatively identical behavior is also found in the gen- 
eral case of non-equal reaction rates k A , ks, and kc (see 
Appendix 0, as the constant of motion (J3J) again guar- 
antees the neutral stability of the deterministic cycles. 



Fluctuations and frequencies around the 
reactive fixed point 



VA 
Vb 



r cos i 
r sin c 



(26) 



Note that in the vicinity of the reactive fixed point, r 
denotes the distance TZ, which is now a random variable: 
TZ = r + o(r). The Fokker-Planck equation (|25|l eventu- 
ally turns into 



d t P(r,<l>,t)=-wod+P(r,<t> i t) 



l r l 



12N 



dl + -d r +dt 



P(r,<f>,t) 



(27) 



The first term on the right-hand-side of this equation de- 
scribes the system's deterministic evolution, being the 
motion on circles around the origin at frequency u>q. 
Stochastic effects enter through the second term, which 
corresponds to isotropic diffusion in two dimensions with 
diffusion constant D = 1/(12AT). Note that it vanishes 
in the limit of infinitely many agents, i.e. when N — > oo. 

If we consider a spherically symmetric probability dis- 
tribution P(r,4>,to) at time to, i.e. independent of the 
angle at to, then this symmetry is conserved by the dy- 
namics according to Eq. I|27[) and one is left with a radial 
distribution function P(r,cf>,t) = P(r,t). The Fokker- 
Planck equation l|27|l thus further simplifies and reads 



d t P(r, t) = D 



rl 



d r + dt 



P(r,t) 



(28) 



This is the diffusion equation in two dimensions with dif- 
fusion constant D = l/(12A r ), expressed in polar coor- 
dinates, for a spherically symmetric probability distri- 
bution. This is the case on which we specifically focus 



In this subsection, we will use the Fokker-Planck equa- 
tion l|28|) to investigate the time-dependence of fluctua- 
tions around the reactive fixed point. In particular, we 
are interested in the time evolution of mean deviation 
from the latter. The average square distance (TZ 2 (t)) will 
be found to grow linearly in time, rescaled by the num- 
ber N of individuals, before saturating. The frequency 
spectrum arising from erratic oscillations around the re- 
active fixed point is of further interest to characterize the 
stochastic dynamics: We will show that the frequency 
|JSJ predicted by the deterministic approach emerges as a 
pole in the power spectrum. 



1. Fluctuations 

Being interested in the vicinity of the reactive fixed 
point, in this subsection we ignore the fact that the 
boundary of the ternary phase space is absorbing. Then, 
the solution to the Fokker-Planck equation (|28|) with 
the initial condition P(r,(f),to = 0) = <J(r), wnere 
the prefactor l/(2?rr) ensures normalization, is simply a 
Gaussian: 



P(r,t) 



1 



■ cxp 



9 



(29) 



4nDt "~ r V 4Dt/ 

This result predicts a broadening of the probability dis- 
tribution in time, as the average square radius increases 
linearly with increasing time t: 



! (t)) = -- 



(30) 



As the time scales linearly with N, increasing the system 
size results in rcscaling the time t and the broadening 
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FIG. 5: (color online) The averaged square radius (7?. 2 ) as 
a function of the rescaled time it = t/N . The blue (or dark 
gray) curve represents a single trajectory, which is seen to 
fluctuate widely. The linear black line indicates the analyt- 
ical prediction 13UH . the red (or light gray) one corresponds 
to sample averages over 10 4 realizations in stochastic simula- 
tions. Hereby, we used a system size of N = 1000. 
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of the probability distribution takes longer. To capture 
these findings, we introduce the rescaled variable u = 
t/N. 

In Fig.|SJ we compare the time evolution of the squared 
radius (lZ 2 (t)) obtained from stochastic simulations with 
the prediction of Eq. (|3U[) and find a good agreement in 
the linear regime around the reactive fixed point, where 
1Z = r + o(r). In fact, for short times, (lZ 2 (t)) displays a 
linear time dependence, with systematic deviation from 
(I30fl at longer times. We understand this as being in part 
due to the linear approximations used to derive i|3U|) , but 
also, and more importantly, to the fact that so far we 
ignored the absorbing character of the boundary. We 
conclude that the latter invalidates the Gaussian prob- 
ability distribution for longer times. This issue, which 
requires a specific analysis, is the scope of Section llll Dl 
where a proper treatment is devised. 

From the finding (r 2 (t)) = 2u/3 together with the 
spherical symmetry of (|29() . resulting in (r(t)<p(t)) = 
((j> 2 {t)} = 0, we readily obtain the variances of the densi- 
ties a and b: 

(x\{t)) = (4(t)) = 

(x A (t)x B (t)) = --u . (31) 

According to these results, the average square deviations 
of the densities from the reactive fixed point grow linearly 
in the rescaled time u, thus exhibiting the same behav- 
ior as we already found for the average squared radius 
(|3Tj)l . In Fig. El these findings are compared to stochastic 
simulations for small times it, where the linear growth is 
indeed recovered. 

We may consider fluctuations of the densities a, b, and 
c not only around the reactive stationary state, as ob- 



FIG. 6: (color online) Time evolution of the variances of the 
densities a and 6 when starting at the reactive fixed point. 
The blue (or dark gray) curves correspond to a single real- 
ization, while the red (or light gray) ones denote averages 
over 1000 samples. Our results are obtained from stochastic 
simulations with a system size of N = 1000. The black line 
indicates the analytical predictions. 



taincd in l|31|) . but also as the deviations from the de- 
terministic cycles. We consider the latter in the linear 
approximation around the fixed point, given by Eq. I|ll|) . 
The fluctuations around them are again described by 
the Fokker-Planck equation l|17(l with a and B given by 
Eqs. (|19f) . but now xa, xb, and xc are the deviations 
from the deterministic cycle: a(t) = ao(t) + XA(t), b(t) = 
b {t) + x B {t), c(t) = co(t) - x A (t) - x B (t), where a , 6 , 
and Co obey Eqs. J5J and characterize the deterministic 
cycles given by Eq. I|ll|) . Again performing van Kam- 
pen's linear noise approximation, we obtain a Fokker- 
Planck equation of the type l|23|) . now with the matrices 

(2a + 2b -l 2a Q \ 

\ -2b Q -2a -2b + l) 

K _ 1 fao{l-a ) -a b \ 

n{ -a b Q bo(l-b )J ■ W 

Note that the entries of the above matrices now depend 
on time t via ao(t) and bo(t). The Fokker-Planck equa- 
tion 123fl yields equations for the time evolution of the 
fluctuations: 

d t (xiXj) = Aik(x k Xj) +Ajk(xkXi) +Bij . (33) 

They may be solved numerically for (x A (t)}, 
(xA{t)xB(t)}, and (xB(t) 2 }, yielding growing oscil- 
lations. In Fig.[7|wc compare these findings to stochastic 




P(w) 



FIG. 7: (color online) Variances of the densities a and b when 
starting at a cycle away from the fixed point. The black lines 
denote our analytical results, while the red (or gray) ones 
are obtained by stochastic simulations as averages over 1 000 
realizations. They are seen to agree for small times, while 
for larger ones the stochasticity of the system induces erratic 
oscillations. These results were obtained when starting from 
a (to = 0) = 0.37, b(t = 0) = 1/3. 



simulations. Wc observe a good agreement for small 
rescaled times u = t/N, while the oscillations of the 
stochastic results become more irregular at longer times. 



2. Frequencies 

Recently, it has been shown that the frequency pre- 
dicted by the deterministic rate equations appears in 
the stochastic system due to a "resonance mechanism" 
|28j . Internal noise present in the system covers all fre- 
quencies and induces excitations; the largest occurring 
for the "resonant" frequency predicted by the rate equa- 
tions. Here, following the same lines, wc address the is- 
sue of the characteristic frequency in the stochastic cyclic 
Lotka-Volterra model. 

In the vicinity of the reactive fixed point, the deter- 
ministic rate equations predict density oscillations 
with frequency cj given in Eq. JSJ. For the stochastic 
model, we now show that a spectrum of frequencies cen- 
tered around this value luq arises. The most convenient 
way to compute this power spectrum P(u>) = (|x(w)| 2 ) 
from the Fokker-Planck equation (|23[) is through the set 
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FIG. 8: The power spectrum for densities in the vicinity of 
the reactive fixed point. The dashed line indicates our analyt- 
ical result, which has a pole at ujo- It agrees with stochastic 
simulations (solid). The inset shows the erratic oscillations 
of the density of one of the species for one realization. The 
system size considered is N = 10000. 



of equivalent Langevin equations [2G 

X{ AijXj ~f~ ^1 



(34) 



with the white noise covariance matrix B: (£,i(t)£,j(t')) — 
BijS(t — t'). From the Fourier transform of Eq. it 
follows that 



P{w) = (|x( W )| 2 ) 
4 1 + 



= Tr[(A- 



3iV (1 - 3w 2 ) 2 



(35) 



The power spectrum has a pole at the characteristic fre- 
quency already predicted from the rate equations J5J), 
luq = 1/V3. For increasing system size N, the power 
spectrum displays a sharper alignment with this value. 
In Fig. [S] we compare our results to stochastic simula- 
tions and find an excellent agreement, except for the pole, 
where the power spectrum in finite systems obviously has 
a finite value. These results were obtained in the vicinity 
of the reactive fixed point, where the linear analysis l|ll|) 
applies. As already found in the deterministic description 
(see Fig. 0}, when departing from the center fixed point, 
nonlincaritics will alter the characteristic frequency. 



D. The extinction probability 

So far, within the stochastic formulation, the fluctua- 
tions around the reactive steady state were found to fol- 
low a Gaussian distribution, linearly broadening in time. 
However, when approaching the absorbing boundary, the 
latter alters this behavior, see Fig. [5J In the following, 
we will incorporate this effect in our quantitative descrip- 
tion. It plays an essential role when discussing the ex- 
tinction probability, which is the scope of this subsection. 
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The probability that one or more species die out in 
the course of time is of special interest within population 
dynamics from a biological viewpoint. When considering 
meta-populations formed of local patches, such questions 
were e.g. raised in Ref. [30l |. Here, we consider the ex- 
tinction probability P ex t(i) that, starting at the reactive 
fixed point, after time t two of the three species have 
died out. In our formulation, this corresponds to the 
probability that after time t the system has reached the 
absorbing boundary of the ternary phase space, depicted 
in Fig. |21 Considering states far from the reactive fixed 
point, we now have to take into account the absorbing 
nature of the boundary. This feature is incorporated in 
our approach by discarding the states having reached the 
boundary, so that a vanishing density of states occurs 
there. As the normalization is lost, we do no longer deal 
with probability distributions: In fact, discarding states 
at the boundary implies a time decay of the integrated 
density of states [which, for commodity, we still refer to 



as P[r,t)\. Ignoring the nonlincarities, the problem now 
takes the form of solving the Fokker-Planck equation H28|) 
for an initial condition P(r,to = 0) = -^:S(r), an d with 
the requirement that P(r, t) has to vanish at the bound- 
ary. Hereby, the triangular shaped absorbing boundary is 
regarded as the outermost (degenerate) cycle of the deter- 
ministic solutions. As the linearization in the y variables 
around the reactive fixed point maps the cycles onto cir- 
cles (see above), the triangular boundary is mapped to 
a sphere as well. Although the linearization scheme on 
which these mappings rely is inaccurate in the vicinity of 
the boundary, it is possible to incorporate nonlinear ef- 
fects in a simple and pragmatic manner, as shown below. 
However, first let us consider the linearized problem, be- 
ing a first-passage to a sphere of radius R, which is, e.g., 
treated in Chapter 6 of Ref. [3l| . The solution is known 
to be a combination of modified Bessel functions of the 
first and second kind. Actually, the Laplace transform of 
the density of states reads 



L. T.{P(r,t)} 



dte 



l P{r,t) 



1 I (R^/D)K (r^/7/D)-I Q (r^7/D)K a (R^/D) 
2tt£> I (Ryfi/D) 



(36) 



where D = l/(12iV) is the diffusion constant; Iq and 
Kq denote the Bessel function of the first, resp. sec- 
ond, kind and of order zero; and R is the radius of the 
absorbing sphere. It is normalized at the initial time 
to = 0; however, for later times t, the total number of 
states J Q R dr J Q 2 ^ dcj> rP(r,t) will decay in time, as states 
are absorbed at the boundary. 

Equipped with this result, we are now in a position 
to calculate the extinction probability P ex t(i). It can be 
found by considering the probability current J(t) at the 
absorbing boundary |3lj . namely: 



J(r) 



D I 

i Sphere 



dP(r, t) 



dr 



r=R 



whose Laplace transform is given by 



L.T.({J(r)}) 



1 



I (R^7/D) 



(37) 



(38) 



The extinction probability at time t is obtained from J(t) 
by integrating over time r until time t |3lj : P ox t(£) = 
J drJ(r). Therefore, the Laplace transform of P ox t{t) 

reads L.T.{P cxt (i)} = l/(sI (R^s/D)). Again, we no- 
tice that we can write this equation in a form that de- 
pends on the time t only via the rescaled time u = t/N: 



L.T.{P ext ( M )} 



sI {R^7jV2) 



(39) 



Hence, for different system sizes N, one obtains the same 
extinction probability P ex t (u) provided one considers the 
same value for the scaling variable t/N = constant (see 
Fig. EJ|. 

We cannot solve the inverse Laplace transform appear- 
ing in Eq. (|39|) . but Iq{z) might be expanded according 
to 113 as 



Io(z) 



(I!) 2 (2!) 2 



±z 2 

4^ 



2\3 



(3!) 5 



(40) 



Considering the first three terms in this expansion, i.e 
contributions up to z 4 , yields the approximate result 



Pext(t) = 1 - 



2 

1 + r 



t 



3R 2 N 



exp 



t 



3R 2 N 



(41) 



Numerically, we have included higher terms. The results 
for contributions up to z 10 in l|40|) are shown in Fig. [5] 
An estimate of the effective distance R to the absorb- 
ing boundary is determined by plugging either a = 0, or 
b = 0, or c = into Eq. 0, which yields R = 1/3. 
From Fig. |U| we observe the extinction probability to 
be overestimated by (|41|l with R = i?o- This stems 
from nonlinearities altering the analysis having led to 
l)4ip. However, adjusting R on physical grounds, we 
are able to capture these effects. Another estimate of 
P GXt is obtained by considering the expression Eq. ijTi|> . 
which arises from a linear analysis. As the extinction of 
two species requires xa = xb = 1/3, in the expression 
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FIG. 9: (color online) The extinction probability when start- 
ing at the reactive fixed point, depending on the rescaled time 
u — t/N. Stochastic simulations for different system sizes 
(N = 100: triangles; N = 200: boxes; N = 500: circles) are 
compared to the analytical prediction 1391 . The left (blue or 
dark gray) is obtained from Rq = 1/3, the right one (red or 
light gray) from Ri = l/y/3, and the middle (black) corre- 
sponds to the average of both: R2 = ^(Ro + Ri). 



Eq. I|14|) this leads to the estimate R\ = 1/V% for R. 
The comparison with Fig. |U| shows that (|41() together 
with R = Ri is a lower bound of P cxt . A simple at- 
tempt to interpolate between the above estimates is to 
consider the mean value of these radii, i?2 = ^(Ro + Ri), 
which happens to yield an excellent agreement with re- 
sults from stochastic simulations, see Fig. |UJ For the 
latter, we have considered systems of N = 100, 200, and 
N = 500 individuals. Rcscaling the time according to 
t/N, they are seen to collapse on a universal curve. This 
is well described by Eq. 1)41(1 with R2 as the radius of the 
absorbing boundary. 

To conclude, we consider the mean time £ a b s that it 
takes until one species becomes extinct. From j^l] we 
find for the rescaled mean absorption time 

u a bs = t ahs /N = 3i? 2 , (42) 

which for our best estimate R — R2 yields a value of 
w abs ~ 0.62. 



IV. CONCLUSION 

Motivated by recent in vitro and in vivo experiments 
aimed at identifying mechanisms responsible for biodi- 
versity in populations of Escherichia coli [Tol Hlj . we 
have considered the stochastic version of the 'rock-paper- 
scissors', or three-species cyclic Lotka-Volterra, system 
within an urn model formulation. This approach allowed 
us to quantitatively study the effect of finite-size fluctu- 
ations in a system with a large, yet finite, number ./V of 
agents. 



While the classical rate equations of the cyclic Lotka- 
Volterra model predict the existence of one (neutrally 
stable) center fixed point, associated with the coexis- 
tence of all the species, this picture is dramatically in- 
validated by the fluctuations which unavoidably appear 
in a finite system. The latter were taken into account by 
a Fokkcr-Planck equation derived from the underlying 
master equation through a Van Kampen expansion. 

Within this scheme, we were able to show that the vari- 
ances of the densities of individuals grow in time (first 
linearly) until extinction of two of the species occurs. 
In this context, we have investigated the probability for 
such extinction to occur at a given time t. As a main 
result of this work, we have shown that this extinction 
probability is a function of the scaling variable t/N. Ex- 
ploiting polar symmetries displayed by the deterministic 
trajectories in the phase portrait and using a mapping 
onto a classical first-passage problem, we were able to 
provide analytic estimates (upper and lower bounds) of 
the extinction probability, which have been successfully 
compared to numerical computation. 

From our results, it turns out that the classical rate 
equation predictions apply to the urn model with a finite 
number of agents only for short enough time, i.e. in the 
regime t -C N. As time increases, the probability of ex- 
tinction grows, asymptotically reaching 1 for t 3> N, so 
that, for finite N, fluctuations are always responsible for 
extinction and thus dramatically jeopardize the possibil- 
ity of coexistence and biodiversity. Interestingly, these 
findings are in qualitative agreement with those (both 
experimental and numerical) reported in Ref . [Tol ] , where 
it was found that in a well mixed environment (as in the 
urn model considered here) two species get extinct. 

While this work has specifically focused on the stochas- 
tic cyclic Lotka-Volterra model, the addressed issues are 
generic. Indeed, we think that our results and techni- 
cal approach, here illustrated by considering the case 
of a paradigmatic model, might actually shed further 
light on the role of fluctuations and the validity of the 
rate equations in a whole class of stochastic systems. In 
fact, while one might believe that fluctuations in an urn 
model should always vanish in the thermodynamic limit, 
we have shown that this issue should be dealt with due 
care: This is true for systems where the rate equations 
predict the existence of an asymptotically stable fixed 
point, which is always reached by the stochastic dynam- 
ics (23,133. In contrast, in systems where the determin- 
istic (rate equation) description predicts the existence of 
(neutrally stable) center fixed points, such as the cyclic 
Lotka-Volterra model, fluctuations have dramatic conse- 
quences and hinder biodiversity by being responsible (at 
long, yet finite, time) for extinction of species [Tol |. In 
this case, instead of a deterministic oscillatory behav- 
ior around the linearly (neutrally) stable fixed point, the 
stochastic dynamics always drives the system toward one 
of its absorbing states. Thus, the absorbing fixed points, 
predicted to be linearly unstable within the rate equation 
theory, actually turn out to be the only stable fixed points 
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available at long time. 



V. ACKNOWLEDGMENTS 

We would like to thank U. C. Tauber and P. L. 
Krapivsky for helpful discussions, as well as A. Traulsen 
and J. C. Clausscn for having made manuscript Q 
available to us prior to publication. M.M. gratefully 
acknowledges the support of the German Alexander 
von Humboldt Foundation through the Fellowship No. 
|iv-scz/1119205|STP. 



APPENDIX A: UNEQUAL REACTION RATES 
RECONSIDERED 

In Section lllll when considering the stochastic ap- 
proach to the cyclic Lotka-Volterra model, for the sake 
of clarity we have specifically turned to the situation of 
equal reaction rates, kA = k B = kc = F In this Ap- 
pendix, we want to provide some details on the general 
case with unequal rates. While the mathematical treat- 
ment becomes more involved, we will argue that the qual- 
itative general situation still follows along the same lines 
as the (simpler) case that we have discussed in detail in 

sec.rrm 

The derivation of the Fokker-Planck equation 12- il) is 
straightforward, following the lines of Subsection IIII Bl 
The matrix A remains unchanged and is given in Eq. |7|>. 
for B we now obtain: 



being no longer proportional to the unit matrix as in the 
case of equal reaction rates. We can do slightly better by 

using an additional rotation, z = ( cos . e o sm V\y, with 

° V — sm V cos 9 J J ' 

rotation angle 



tan(26>) = 



3(fc c - k A )\Jk A kBkc{kA + k B + k c ) 
k\{k c - k B ) + k 2 c {k A - k B ) 



(A4) 

This variable transformation leaves A invariant, but 
brings B to diagonal form, with unequal diagonal el- 
ements. The stochastic effects thus correspond to 
anisotropic diffusion. However, for large system size N, 
the effects of the anisotropy are washed out: The sys- 
tem's motion on the deterministic cycles, described by 
A, occurs on a much faster timescale then the anisotropic 
diffusion, resulting in an averaging over the different di- 
rections. 

To calculate the time evolution of the average deviation 
from the reactive fixed point, dt (1Z 2 (t)), we start from the 
the fluctuations in y, which fulfill the equations 



dtiUiVj) = A lk {vkVj) + AjkiVkVi) + Br 



Using Eq. (|15f) . we obtain 



(A5) 



d t {K 2 (t)) =Baa + Bbb 



(A6) 



B 



1 



N k A + k B + k c 



2 -1 
-1 2 



(Al) 



The corresponding Fokker-Planck equation reads 



9 t P(x, t) = -dilAijXjPix, t)] + ^BijdidjPix, f) . 

(A2) 

Again, we aim at benefitting from the cyclic structure of 
the deterministic solutions, and perform a variable trans- 
formation to y = <Sx, with S given in Eq. @ ■ As already 

found in Subsection 111 Bl A turns into A = f_^ q)> 

such that in the y variables the deterministic solutions 
correspond to circles around the origin. For the stochas- 
tic part, entering via B, a technical difficulty arises. It is 
transformed into 



B=l l 



2 N k A + k B + k c 



' k~ A +k A k c +k c _ 9 k A -k c 
k^k% W 2k A k c 

1 



i—&c , , 



UJ 



(A3) 



The dependence on the deterministic part has dropped 
out, and the solution to the above equation with initial 
condition (lZ 2 (t = 0)) = is a linear increase in the 
rescaled time u = t/N: 



k A k 



A^c 



^o 2 +l 



t 

2k A + k B + k c [ fep| ~ U ° +1 \n ' 

(A7) 

valid around the reactive fixed point. As in the case 
of equal reaction rates, we have a linear dependence 
(lZ 2 (t)) ~ t/N, corresponding to a two-dimensional ran- 
dom walk. 

As a conclusion of the above discussion, the general 
case of unequal reaction rates qualitatively reproduces 
the behavior of the before discussed simplest situation of 
equal rates (confirmed by stochastic simulations). The 
latter turns out to already provide a comprehensive un- 
derstanding of the system. 
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